# Cargar tmap
library(tmap)
# Colección de paquetes de Tidyverse
library(tidyverse)
# Estilos para ggplot2
library(ggthemes)
# Paletas de colores de RColorBrewer
library(RColorBrewer)
# Paletas de colores de viridis
library(viridisLite)
# Gráficos interactivos
library(plotly)
# Manejo de datos vectoriales
library(sf)
# Manejo de datos raster
library(terra)
# Manejo de datos raster
library(raster)
# Mapas interactivos
library(leaflet)
# Acceso a datos en GBIF
library(rgbif)
# Datos geoespaciales
library(geodata)
# Modelado de distribución de especies
library(dismo)
# Manejo datos Maxent
library(rJava)
# Cargar la librería de widgets
library(htmlwidgets)Aves Guanacaste
Cargar librerías
Contexto
Los datos de este estudio se obtuvieron de tres fuentes distintas:
eBird: que es una plataforma de ciencia ciudadana que recopila, almacena y analiza datos de aves de personas observadoras en todo el mundo.
WorldClim: proporciona datos climáticos a nivel mundial, estos incluyen variable climáticas como tempertura, precipitación, radiación solar, humedad relativa y otros factores importantes para la modelación climatológica y ecológica.
HydroSHED: que tiene datos acerca del flujo del agua, acumulación de flujo, cuencas hidrológicas, río y lagos y permiten el análisis detallado sobre los sistemas hidrográficos a escala global.
Definimos el área de Guanacaste como área de estudio debido a que es un territorio que cuenta con 2 sitios Ramsar: el Embalse Arenal (8,317 ha) y el humedal del PN Palo Verde los cuales son sitios de reproducción y alimentación para una gran cantidad de especies de aves acuáticas, migratorias y residentes. Estos humedales funcionan como refugio de especies en peligro de extinción y constituye una de las zonas de anidamiento más grandes del país (SINAC,2020 y 2024)
Nuestro objeto de estudio son las aves: Dendrocygna autumnalis, Spatula discors, Jabiru mycteria, Mycteria americana, Eudocimus albus. Estas especies tienen una vida vinculada a los humedales, pues sus poblaciones y ciclos de vida dependen de estos ecosistemas (alimentación, reproducción, anidación, etc). Todas son residentes a excepción Spatula discors que tiene un periodo de migración al norte, entre septiembre y abril, para reproducirse (Barchiesi et al., 2021).
Cargar datos aves y área de estudio
# Lectura de un archivo TXT con registros de presencia de aves en Costa Rica
puntos <- read.delim(
"C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Aves_Palo_Verde/Ebird/Ebird data Guanacaste/ebd_CR-G_200701_202412_unv_smp_relSep-2024/ebd_CR-G_200701_202412_unv_smp_relSep-2024.txt"
)
# Cargar shape Guanacaste
guanacaste <- st_read("C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Monitoreo_Palo_Verde/Guanacaste/delimitshapes/guanacaste.shp",
quiet = TRUE)Filtrar especies seleccionadas
# Filtrar las especies de aves de interés
especies_seleccionadas <- puntos |>
filter(SCIENTIFIC.NAME == "Dendrocygna autumnalis" |
SCIENTIFIC.NAME == "Spatula discors" |
SCIENTIFIC.NAME == "Jabiru mycteria" |
SCIENTIFIC.NAME == "Mycteria americana" |
SCIENTIFIC.NAME == "Eudocimus albus")Definición espacial: sistema referenciado de coordenadas.
# Convertir a objeto espacial usando columnas de coordenadas
especies_seleccionadas <- st_as_sf(
especies_seleccionadas,
coords = c("LONGITUDE", "LATITUDE"), # Nombres exactos de las columnas
remove = FALSE, # conservar las columnas de las coordenadas
crs = 4326 # Sistema de referencia WGS84 (EPSG:4326)
)
# Asignación de un CRS al objeto guanacaste
guanacaste <-
guanacaste |>
st_transform(4326)#Dar formato de fecha a la columna correspondiente
especies_seleccionadas$OBSERVATION.DATE <- as.Date(especies_seleccionadas$OBSERVATION.DATE, format = "%Y-%m-%d")Graficar registros de presencia por año por especie
Podemos observar la cantidad de registros de las especies dependen del interés de las y los observadores, por lo que especies interesantes o difíciles de conseguir presentan muchos registros, como el Mycteria americana, a diferencia de aquellas que son más especies más comunes y por tanto más fáciles de encontrar como Dendrocygna autumnalis
# Gráfico ggplot2
grafico_ggplot2 <-
especies_seleccionadas |>
st_drop_geometry() |>
group_by(year = year(OBSERVATION.DATE), SCIENTIFIC.NAME) |>
summarize(n = n()) |>
ggplot(aes(x = year, y = n, color=SCIENTIFIC.NAME)) +
geom_line() +
geom_point(
aes(
text = paste0(
"Año: ", year, "\n",
"Cantidad de registros: ", n
)
)
) +
ggtitle("Cantidad de registros de presencia por año") +
xlab("Año") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: eBird") +
theme_minimal()
# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |>
config(locale = 'es')#Gráfico de barras anual Se muestra la coantidad de individuos por especie registrados cada año. Se puede observa que los patos son quienes presentan mayor número de individuos (Dendrocygna autumnalis y Spatula discors).
# Agrupo y cuento los individuos por año
df <- especies_seleccionadas |>
group_by(SCIENTIFIC.NAME, year = year(OBSERVATION.DATE)) |>
summarise(
individuos = sum(as.numeric(OBSERVATION.COUNT), na.rm = TRUE)
) # Elimino nans
# Gráfico de barras
graf_barras <- df |>
ggplot(aes(x = as.factor(year), y = individuos, fill = as.factor(SCIENTIFIC.NAME))) +
geom_bar(
stat = "identity",
aes(
text = paste0(
" Individuos por año: ", round(individuos, 2)
)
)
) +
ggtitle("Número de Individuos por Año") +
xlab("Año") +
ylab("Cantidad de Individuos") +
labs(caption = "Fuente: Datos de eBird", fill = "Especies") +
theme_minimal()
# scale_fill_manual(values = colores, name = "Año") # Pongo los colores definidos y nombre a la leyenda
# Gráfico interactivo con Plotly
ggplotly(graf_barras, tooltip = "text") |>
config(locale = 'es')Cargar los archivos de Clima
# Definir el directorio de los datos climáticos
ruta_clima <- "C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Aves_Palo_Verde/Modelos/climate/bios/"
# Listar los archivos .tif en el directorio
archivos_clima <- list.files(ruta_clima, pattern = "\\.tif$", full.names = TRUE)
# Cargar todos los archivos raster
bios <- rast(archivos_clima)
names(bios) [1] "wc2.1_30s_bio_1" "wc2.1_30s_bio_10" "wc2.1_30s_bio_11" "wc2.1_30s_bio_12"
[5] "wc2.1_30s_bio_13" "wc2.1_30s_bio_14" "wc2.1_30s_bio_15" "wc2.1_30s_bio_16"
[9] "wc2.1_30s_bio_17" "wc2.1_30s_bio_18" "wc2.1_30s_bio_19" "wc2.1_30s_bio_2"
[13] "wc2.1_30s_bio_3" "wc2.1_30s_bio_4" "wc2.1_30s_bio_5" "wc2.1_30s_bio_6"
[17] "wc2.1_30s_bio_7" "wc2.1_30s_bio_8" "wc2.1_30s_bio_9" Recortar área de estudio
# Definir la extensión del área de estudio
area_estudio <- ext(
min(especies_seleccionadas$LONGITUDE) - 5,
max(especies_seleccionadas$LONGITUDE) + 5,
min(especies_seleccionadas$LATITUDE) - 5,
max(especies_seleccionadas$LATITUDE) + 5
)
# Recortar las variables bioclimáticas al área de estudio
clima <- crop(bios, area_estudio)Cargar datos de humedales
Nosotras utilizamos la capa de Global Lakes and Wetlands Database: GLWD_v2_delta_area_ha_x10, que nos indica el porcentaje de cobertura de humedales lo cual para éste análisis es adecuado puesto que las aves sus hábitos de vida están relacionados a ellos.
ruta_humedal <- "C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Aves_Palo_Verde/GLWD_v2_delta_combined_classes_tif/GLWD_v2_delta_combined_classes/GLWD_v2_delta_area_ha_x10.tif"
archivo_humedales <- rast(ruta_humedal)
# Recortar y resamplear 'humedales' para que tenga la misma extensión y resolución que 'clima'
humedales <- archivo_humedales |>
crop(ext(clima)) |>
resample(clima)Unir datos de clima y humedales en un mismo ráster
# Unir las capas de clima y de humedales
datos <- c(clima, humedales)Entrenamiento del modelo MaxEnt
Se crea un DataFrame con las coordenadas de Longitud y Latitud
# Crear dataframe con columnas de longitud y latitud
coordenadas_presencia <- data.frame(
decimalLongitude = especies_seleccionadas$LONGITUDE,
decimalLatitude = especies_seleccionadas$LATITUDE
)
# Eliminar coordenadas duplicadas
coordenadas_presencia <- unique(coordenadas_presencia)Crear semilla para garantizar aleatoriedad reproducible
# Establecer una "semilla" para garantizar que la selección aleatoria sea reproducible
set.seed(123)
# Cantidad de registros de presencia
n_presencia <- nrow(coordenadas_presencia)
# Con sample(), se selecciona aleatoriamente una proporción (ej. 0.7)
# de los índices de los datos de presencia para el conjunto de entrenamiento
indices_entrenamiento <- sample(
1:n_presencia,
size = round(0.7 * n_presencia)
)
# Crear el subconjunto de entrenamiento utilizando los índices seleccionados
entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]
# Crear el subconjunto de evaluación con los datos restantes
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]Cambiar formato a ráster y ejecutar el modelo
# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a que es este el que acepta el paquete dismo
datos <- raster::stack(datos)
# Ejecutar el modelo
modelo_maxent <- maxent(x = datos, p = entrenamiento)
# Aplicar el modelo entrenado a las variables climáticas
# para generar un mapa de idoneidad del hábitat
prediccion <- predict(modelo_maxent, datos)Evaluacion del modelo MaxEnt
# terra::extract() extrae los valores del raster de predicción
# en las coordenadas de evaluación
# eval_pres almacena los valores de idoneidad predichos
# en los puntos de evaluación de presencia
eval_pres <- terra::extract(
prediccion,
evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)
# Generar puntos aleatorios dentro del área de estudio definida.
# Estos puntos se asumen como ausencias de la especie.
ausencias <- randomPoints(mask = datos, n = 1000)
# eval_aus almacena los valores de idoneidad predichos
# en los puntos de ausencia
eval_aus <- terra::extract(
prediccion,
ausencias
)
# Generar estadísticas de evaluación del modelo
resultado_evaluacion <- evaluate(p = eval_pres, a = eval_aus)Curva ROC
# Datos para graficar la curva ROC
datos_roc <- data.frame(
FPR = resultado_evaluacion@FPR,
TPR = resultado_evaluacion@TPR,
Umbral = resultado_evaluacion@t
)
# Valor AUC
auc <- resultado_evaluacion@auc
# Gráfico ggplot2
grafico_ggplot2 <-
ggplot(
datos_roc,
aes(
x = FPR,
y = TPR,
u = Umbral
)
) +
geom_line(
color = "blue",
size = 1
) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
labs(title = paste("Curva ROC (AUC =", round(auc, 3), ")"),
x = "Tasa de falsos positivos (FPR)",
y = "Tasa de verdaderos positivos (TPR)") +
theme_minimal()Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.# Gráfico plotly
ggplotly(grafico_ggplot2) |>
config(locale = 'es')#Graficar mapa de distribución continua
# Paleta de colores del modelo
colores_modelo <- colorNumeric(
palette = c("white", "black"),
values(prediccion),
na.color = "transparent"
)
# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
# palette = "inferno",
# palette = "magma",
palette = rev(brewer.pal(11, "RdYlBu")),
values(datos$wc2.1_30s_bio_1),
na.color = "transparent"
)
# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
# palette = "viridis",
# palette = "YlGnBu",
palette = "Blues",
values(datos$wc2.1_30s_bio_12),
na.color = "transparent"
)
# Paleta de colores de precipitación
colores_humedales <- colorNumeric(
# palette = "viridis",
# palette = "YlGnBu",
palette = "YlGnBu",
values(datos$Band_1),
na.color = "transparent"
)
# Crear una paleta de colores categórica para las especies
especies_unicas <- unique(especies_seleccionadas$SCIENTIFIC.NAME)
paleta_especies <- colorFactor(
palette = brewer.pal(length(especies_unicas), "Set1"),
domain = especies_unicas
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
datos$wc2.1_30s_bio_1,
colors = colores_temperatura, # paleta de colores
opacity = 0.6,
group = "Temperatura",
) |>
addRasterImage( # capa raster de precipitación
datos$wc2.1_30s_bio_12,
colors = colores_precipitacion, # paleta de colores
opacity = 0.6,
group = "Precipitación",
) |>
addRasterImage( # capa raster de humedales
datos$Band_1,
colors = colores_humedales, # paleta de colores
opacity = 0.6,
group = "Humedales",
) |>
addRasterImage( # capa raster del modelo de distribución
prediccion,
colors = colores_modelo,
opacity = 0.6,
group = "Modelo de distribución",
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = especies_seleccionadas,
stroke = F,
radius = 3,
fillColor = ~paleta_especies(SCIENTIFIC.NAME),
fillOpacity = 1,
popup = paste(
paste0("<strong>Nombre Científico: </strong>", especies_seleccionadas$SCIENTIFIC.NAME),
paste0("<strong>Cantidad de individuos: </strong>", especies_seleccionadas$OBSERVATION.COUNT),
paste0("<strong>Fecha de avistamiento: </strong>", especies_seleccionadas$OBSERVATION.DATE),
paste0("<strong>Fecha de avistamiento: </strong>", especies_seleccionadas$COMMON.NAME),
paste0("<a href='", especies_seleccionadas$GLOBAL.UNIQUE.IDENTIFIER, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de aves"
) |>
addLegend(
title = "Temperatura",
values = values(datos$wc2.1_30s_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(datos$wc2.1_30s_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Humedales",
values = values(datos$Band_1),
pal = colores_humedales,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Modelo de distribución",
values = values(prediccion),
pal = colores_modelo,
position = "bottomright",
group = "Modelo de distribución"
) |>
addLayersControl(
# control de capas
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c(
"Temperatura",
"Precipitación",
"Humedales",
"Modelo de distribución",
"Registros de aves"
)
) |>
hideGroup("Temperatura") |>
hideGroup("Precipitación") |>
hideGroup("Humedales")Warning in colors(.): Some values were outside the color scale and will be
treated as NA